home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / expint.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  11.6 KB  |  449 lines
  1. /* specfunc/expint.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author: G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_expint.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "chebyshev.h"
  30. #include "cheb_eval.c"
  31.  
  32. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  33.  
  34. /*
  35.  Chebyshev expansions: based on SLATEC e1.f, W. Fullerton
  36.  
  37.  Series for AE11       on the interval -1.00000D-01 to  0.
  38.                     with weighted error   1.76E-17
  39.                      log weighted error  16.75
  40.                    significant figures required  15.70
  41.                     decimal places required  17.55
  42.  
  43.  
  44.  Series for AE12       on the interval -2.50000D-01 to -1.00000D-01
  45.                     with weighted error   5.83E-17
  46.                      log weighted error  16.23
  47.                    significant figures required  15.76
  48.                     decimal places required  16.93
  49.  
  50.  
  51.  Series for E11        on the interval -4.00000D+00 to -1.00000D+00
  52.                     with weighted error   1.08E-18
  53.                      log weighted error  17.97
  54.                    significant figures required  19.02
  55.                     decimal places required  18.61
  56.  
  57.  
  58.  Series for E12        on the interval -1.00000D+00 to  1.00000D+00
  59.                     with weighted error   3.15E-18
  60.                      log weighted error  17.50
  61.             approx significant figures required  15.8
  62.                     decimal places required  18.10
  63.  
  64.  
  65.  Series for AE13       on the interval  2.50000D-01 to  1.00000D+00
  66.                     with weighted error   2.34E-17
  67.                      log weighted error  16.63
  68.                    significant figures required  16.14
  69.                     decimal places required  17.33
  70.  
  71.  
  72.  Series for AE14       on the interval  0.        to  2.50000D-01
  73.                     with weighted error   5.41E-17
  74.                      log weighted error  16.27
  75.                    significant figures required  15.38
  76.                     decimal places required  16.97
  77. */
  78.  
  79. static double AE11_data[39] = {
  80.    0.121503239716065790,
  81.   -0.065088778513550150,
  82.    0.004897651357459670,
  83.   -0.000649237843027216,
  84.    0.000093840434587471,
  85.    0.000000420236380882,
  86.   -0.000008113374735904,
  87.    0.000002804247688663,
  88.    0.000000056487164441,
  89.   -0.000000344809174450,
  90.    0.000000058209273578,
  91.    0.000000038711426349,
  92.   -0.000000012453235014,
  93.   -0.000000005118504888,
  94.    0.000000002148771527,
  95.    0.000000000868459898,
  96.   -0.000000000343650105,
  97.   -0.000000000179796603,
  98.    0.000000000047442060,
  99.    0.000000000040423282,
  100.   -0.000000000003543928,
  101.   -0.000000000008853444,
  102.   -0.000000000000960151,
  103.    0.000000000001692921,
  104.    0.000000000000607990,
  105.   -0.000000000000224338,
  106.   -0.000000000000200327,
  107.   -0.000000000000006246,
  108.    0.000000000000045571,
  109.    0.000000000000016383,
  110.   -0.000000000000005561,
  111.   -0.000000000000006074,
  112.   -0.000000000000000862,
  113.    0.000000000000001223,
  114.    0.000000000000000716,
  115.   -0.000000000000000024,
  116.   -0.000000000000000201,
  117.   -0.000000000000000082,
  118.    0.000000000000000017
  119. };
  120. static cheb_series AE11_cs = {
  121.   AE11_data,
  122.   38,
  123.   -1, 1,
  124.   20
  125. };
  126.  
  127. static double AE12_data[25] = {
  128.    0.582417495134726740,
  129.   -0.158348850905782750,
  130.   -0.006764275590323141,
  131.    0.005125843950185725,
  132.    0.000435232492169391,
  133.   -0.000143613366305483,
  134.   -0.000041801320556301,
  135.   -0.000002713395758640,
  136.    0.000001151381913647,
  137.    0.000000420650022012,
  138.    0.000000066581901391,
  139.    0.000000000662143777,
  140.   -0.000000002844104870,
  141.   -0.000000000940724197,
  142.   -0.000000000177476602,
  143.   -0.000000000015830222,
  144.    0.000000000002905732,
  145.    0.000000000001769356,
  146.    0.000000000000492735,
  147.    0.000000000000093709,
  148.    0.000000000000010707,
  149.   -0.000000000000000537,
  150.   -0.000000000000000716,
  151.   -0.000000000000000244,
  152.   -0.000000000000000058
  153. };
  154. static cheb_series AE12_cs = {
  155.   AE12_data,
  156.   24,
  157.   -1, 1,
  158.   15
  159. };
  160.  
  161. static double E11_data[19] = {
  162.   -16.11346165557149402600,
  163.     7.79407277874268027690,
  164.    -1.95540581886314195070,
  165.     0.37337293866277945612,
  166.    -0.05692503191092901938,
  167.     0.00721107776966009185,
  168.    -0.00078104901449841593,
  169.     0.00007388093356262168,
  170.    -0.00000620286187580820,
  171.     0.00000046816002303176,
  172.    -0.00000003209288853329,
  173.     0.00000000201519974874,
  174.    -0.00000000011673686816,
  175.     0.00000000000627627066,
  176.    -0.00000000000031481541,
  177.     0.00000000000001479904,
  178.    -0.00000000000000065457,
  179.     0.00000000000000002733,
  180.    -0.00000000000000000108
  181. };
  182. static cheb_series E11_cs = {
  183.   E11_data,
  184.   18,
  185.   -1, 1,
  186.   13
  187. };
  188.  
  189. static double E12_data[16] = {
  190.   -0.03739021479220279500,
  191.    0.04272398606220957700,
  192.   -0.13031820798497005440,
  193.    0.01441912402469889073,
  194.   -0.00134617078051068022,
  195.    0.00010731029253063780,
  196.   -0.00000742999951611943,
  197.    0.00000045377325690753,
  198.   -0.00000002476417211390,
  199.    0.00000000122076581374,
  200.   -0.00000000005485141480,
  201.    0.00000000000226362142,
  202.   -0.00000000000008635897,
  203.    0.00000000000000306291,
  204.   -0.00000000000000010148,
  205.    0.00000000000000000315
  206. };
  207. static cheb_series E12_cs = {
  208.   E12_data,
  209.   15,
  210.   -1, 1,
  211.   10
  212. };
  213.  
  214. static double AE13_data[25] = {
  215.   -0.605773246640603460,
  216.   -0.112535243483660900,
  217.    0.013432266247902779,
  218.   -0.001926845187381145,
  219.    0.000309118337720603,
  220.   -0.000053564132129618,
  221.    0.000009827812880247,
  222.   -0.000001885368984916,
  223.    0.000000374943193568,
  224.   -0.000000076823455870,
  225.    0.000000016143270567,
  226.   -0.000000003466802211,
  227.    0.000000000758754209,
  228.   -0.000000000168864333,
  229.    0.000000000038145706,
  230.   -0.000000000008733026,
  231.    0.000000000002023672,
  232.   -0.000000000000474132,
  233.    0.000000000000112211,
  234.   -0.000000000000026804,
  235.    0.000000000000006457,
  236.   -0.000000000000001568,
  237.    0.000000000000000383,
  238.   -0.000000000000000094,
  239.    0.000000000000000023
  240. };
  241. static cheb_series AE13_cs = {
  242.   AE13_data,
  243.   24,
  244.   -1, 1,
  245.   15
  246. };
  247.  
  248. static double AE14_data[26] = {
  249.   -0.18929180007530170,
  250.   -0.08648117855259871,
  251.    0.00722410154374659,
  252.   -0.00080975594575573,
  253.    0.00010999134432661,
  254.   -0.00001717332998937,
  255.    0.00000298562751447,
  256.   -0.00000056596491457,
  257.    0.00000011526808397,
  258.   -0.00000002495030440,
  259.    0.00000000569232420,
  260.   -0.00000000135995766,
  261.    0.00000000033846628,
  262.   -0.00000000008737853,
  263.    0.00000000002331588,
  264.   -0.00000000000641148,
  265.    0.00000000000181224,
  266.   -0.00000000000052538,
  267.    0.00000000000015592,
  268.   -0.00000000000004729,
  269.    0.00000000000001463,
  270.   -0.00000000000000461,
  271.    0.00000000000000148,
  272.   -0.00000000000000048,
  273.    0.00000000000000016,
  274.   -0.00000000000000005
  275. };
  276. static cheb_series AE14_cs = {
  277.   AE14_data,
  278.   25,
  279.   -1, 1,
  280.   13
  281. };
  282.  
  283.  
  284. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  285.  
  286.  
  287. int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
  288. {
  289.   const double xmaxt = -GSL_LOG_DBL_MIN;      /* XMAXT = -LOG (R1MACH(1)) */
  290.   const double xmax  = xmaxt - log(xmaxt);    /* XMAX = XMAXT - LOG(XMAXT) */
  291.  
  292.   /* CHECK_POINTER(result) */
  293.  
  294.   if(x < -xmax) {
  295.     OVERFLOW_ERROR(result);
  296.   }
  297.   else if(x <= -10.0) {
  298.     const double s = exp(-x)/x;
  299.     gsl_sf_result result_c;
  300.     cheb_eval_e(&AE11_cs, 20.0/x+1.0, &result_c);
  301.     result->val  = s * (1.0 + result_c.val);
  302.     result->err  = s * result_c.err;
  303.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) * fabs(result->val);
  304.     return GSL_SUCCESS;
  305.   }
  306.   else if(x <= -4.0) {
  307.     const double s = exp(-x)/x;
  308.     gsl_sf_result result_c;
  309.     cheb_eval_e(&AE12_cs, (40.0/x+7.0)/3.0, &result_c);
  310.     result->val  = s * (1.0 + result_c.val);
  311.     result->err  = s * result_c.err;
  312.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  313.     return GSL_SUCCESS;
  314.   }
  315.   else if(x <= -1.0) {
  316.     const double ln_term = -log(fabs(x));
  317.     gsl_sf_result result_c;
  318.     cheb_eval_e(&E11_cs, (2.0*x+5.0)/3.0, &result_c);
  319.     result->val  = ln_term + result_c.val;
  320.     result->err  = result_c.err + GSL_DBL_EPSILON * fabs(ln_term);
  321.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  322.     return GSL_SUCCESS;
  323.   }
  324.   else if(x == 0.0) {
  325.     DOMAIN_ERROR(result);
  326.   }
  327.   else if(x <= 1.0) {
  328.     const double ln_term = -log(fabs(x));
  329.     gsl_sf_result result_c;
  330.     cheb_eval_e(&E12_cs, x, &result_c);
  331.     result->val  = ln_term - 0.6875 + x + result_c.val;
  332.     result->err  = result_c.err + GSL_DBL_EPSILON * fabs(ln_term);
  333.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  334.     return GSL_SUCCESS;
  335.   }
  336.   else if(x <= 4.0) {
  337.     const double s = exp(-x)/x;
  338.     gsl_sf_result result_c;
  339.     cheb_eval_e(&AE13_cs, (8.0/x-5.0)/3.0, &result_c);
  340.     result->val  = s * (1.0 + result_c.val);
  341.     result->err  = s * result_c.err;
  342.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  343.     return GSL_SUCCESS;
  344.   }
  345.   else if(x <= xmax) {
  346.     const double s = exp(-x)/x;
  347.     gsl_sf_result result_c;
  348.     cheb_eval_e(&AE14_cs, 8.0/x-1.0, &result_c);
  349.     result->val  = s * (1.0 +  result_c.val);
  350.     result->err  = s * (GSL_DBL_EPSILON + result_c.err);
  351.     result->err += 2.0 * (x + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
  352.     return GSL_SUCCESS;
  353.   }
  354.   else {
  355.     UNDERFLOW_ERROR(result);
  356.   }
  357. }
  358.  
  359.  
  360. int gsl_sf_expint_E2_e(const double x, gsl_sf_result * result)
  361. {
  362.   const double xmaxt = -GSL_LOG_DBL_MIN;
  363.   const double xmax  = xmaxt - log(xmaxt);
  364.  
  365.   /* CHECK_POINTER(result) */
  366.  
  367.   if(x < -xmax) {
  368.     OVERFLOW_ERROR(result);
  369.   }
  370.   else if(x < 100.0) {
  371.     const double ex = exp(-x);
  372.     gsl_sf_result result_E1;
  373.     int stat_E1 = gsl_sf_expint_E1_e(x, &result_E1);
  374.     result->val  = ex - x*result_E1.val;
  375.     result->err  = fabs(x) * (GSL_DBL_EPSILON*ex + result_E1.err);
  376.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  377.     return stat_E1;
  378.   }
  379.   else if(x < xmax) {
  380.     const double c1  = -2.0;
  381.     const double c2  =  6.0;
  382.     const double c3  = -24.0;
  383.     const double c4  =  120.0;
  384.     const double c5  = -720.0;
  385.     const double c6  =  5040.0;
  386.     const double c7  = -40320.0;
  387.     const double c8  =  362880.0;
  388.     const double c9  = -3628800.0;
  389.     const double c10 =  39916800.0;
  390.     const double c11 = -479001600.0;
  391.     const double c12 =  6227020800.0;
  392.     const double c13 = -87178291200.0;
  393.     const double y = 1.0/x;
  394.     const double sum6 = c6+y*(c7+y*(c8+y*(c9+y*(c10+y*(c11+y*(c12+y*c13))))));
  395.     const double sum  = y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*sum6)))));
  396.     result->val = exp(-x) * (1.0 + sum)/x;
  397.     result->err = 2.0 * (x + 1.0) * GSL_DBL_EPSILON * result->val;
  398.     return GSL_SUCCESS;
  399.   }
  400.   else {
  401.     UNDERFLOW_ERROR(result);
  402.   }
  403. }
  404.  
  405.  
  406. int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result)
  407. {
  408.   /* CHECK_POINTER(result) */
  409.  
  410.   {
  411.     int status = gsl_sf_expint_E1_e(-x, result);
  412.     result->val = -result->val;
  413.     return status;
  414.   }
  415. }
  416.  
  417. #if 0
  418. static double recurse_En(int n, double x, double E1)
  419. {
  420.   int i;
  421.   double En = E1;
  422.   double ex = exp(-x);
  423.   for(i=2; i<=n; i++) {
  424.     En = 1./(i-1) * (ex - x * En);
  425.   }
  426.   return En;
  427. }
  428. #endif
  429.  
  430.  
  431. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  432.  
  433. #include "eval.h"
  434.  
  435. double gsl_sf_expint_E1(const double x)
  436. {
  437.   EVAL_RESULT(gsl_sf_expint_E1_e(x, &result));
  438. }
  439.  
  440. double gsl_sf_expint_E2(const double x)
  441. {
  442.   EVAL_RESULT(gsl_sf_expint_E2_e(x, &result));
  443. }
  444.  
  445. double gsl_sf_expint_Ei(const double x)
  446. {
  447.   EVAL_RESULT(gsl_sf_expint_Ei_e(x, &result));
  448. }
  449.